

a=1;
t0=0;
T=1.0;
f=@(x,t) a*x;
ft=@(x,t) 0;
fx=@(x,t) a;

sol=@(t) exp(a*(t-t0));
x0=1.0; 

N=20;
M=8;
er=zeros(M,5);
for k=1:M,
 [x,t]=exEuler19(f,t0,T,x0,N);
 x1=exp(a*(T-t0)/N);
 [xm]=midpoint19(f,t0,T,x0,x1,N);
 [xt]=taylor19(f,ft,fx,x0,t0,T,N);
 [xa,t]=AdamsBashforth19(f,t0,T,x0,x1,N);
 [xH]=Heun19(f,t0,T,x0,N);
 N*=2;
 %er(k,1)=sol(T)-x(end);
 er(k,1)=norm(sol(t)-x,'inf');
 %er(k,2)=sol(T)-xm(end);
 er(k,2)=norm(sol(t)-xm,'inf');
 %er(k,3)=sol(T)-xt(end);
 er(k,3)=norm(sol(t)-xt,'inf');
  %er(k,4)=sol(T)-xa(end);
 er(k,4)=norm(sol(t)-xa,'inf');
  %er(k,5)=sol(T)-xH(end);
 er(k,5)=norm(sol(t)-xH,'inf');
end 
er
er(1:M-1,:)./er(2:M,:)


%dx/dt=-x

a=-1;
t0=0;
T=1.0;
f=@(x,t) a*x;
ft=@(x,t) 0;
fx=@(x,t) a;

sol=@(t) exp(a*(t-t0));
x0=1.0; 

N=20;
M=8;
er=zeros(M,5);
for k=1:M,
 [x,t]=exEuler19(f,t0,T,x0,N);
 x1=exp(a*(T-t0)/N);
 [xm]=midpoint19(f,t0,T,x0,x1,N);
 [xt]=taylor19(f,ft,fx,x0,t0,T,N);
 [xa,t]=AdamsBashforth19(f,t0,T,x0,x1,N);
 [xH]=Heun19(f,t0,T,x0,N);
 N*=2;
 %er(k,1)=sol(T)-x(end);
 er(k,1)=norm(sol(t)-x,'inf');
 %er(k,2)=sol(T)-xm(end);
 er(k,2)=norm(sol(t)-xm,'inf');
 %er(k,3)=sol(T)-xt(end);
 er(k,3)=norm(sol(t)-xt,'inf');
  %er(k,4)=sol(T)-xa(end);
 er(k,4)=norm(sol(t)-xa,'inf');
  %er(k,5)=sol(T)-xH(end);
 er(k,5)=norm(sol(t)-xH,'inf');
end 
er
er(1:M-1,:)./er(2:M,:)


%tests for penduulum equartion (linear)
%y''=-y y(0)=1;y'(0)=0;
%y'(1)=y(2) y(1)(t=0)=1
%y'(2)=-y(1) y(2)(t=0)=0
a=[0,1;-1,0];
%y'=a*y
t0=0;
T=1.0;
f=@(x,t) a*x;
ft=@(x,t) zeros(size(x));
fx=@(x,t) a;

sol=@(t) [cos(t);-sin(t)]
x0=[1;0] 

N=20;
M=8;
er=zeros(M,5);
for k=1:M,
 [x,t]=exEuler19(f,t0,T,x0,N);
 h=(T-t0)/N;x1=sol(h);
 [xm]=midpoint19(f,t0,T,x0,x1,N);
 [xt]=taylor19(f,ft,fx,x0,t0,T,N);
 [xa,t]=AdamsBashforth19(f,t0,T,x0,x1,N);
 [xH]=Heun19(f,t0,T,x0,N);
 N*=2;
% er(k,1)=norm(sol(T)-x(end,:)');
 er(k,1)=norm(sol(t')'-x,"inf");
% er(k,2)=norm(sol(T)-xm(end,:)');
 er(k,2)=norm(sol(t')'-xm,"inf");
%er(k,3)=norm(sol(T)-xt(end,:)');
 er(k,3)=norm(sol(t')'-xt,"inf");
 %er(k,4)=norm(sol(T)-xa(end,:)');
 er(k,4)=norm(sol(t')'-xa,"inf");
 %er(k,5)=norm(sol(T)-xH(end,:)');
 er(k,5)=norm(sol(t')'-xH,"inf");
end 
er
er(1:M-1,:)./er(2:M,:)



